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I.  Abstract 

This  paper  focuses  on  the  time  integration  of  the  nonlinear  EOM  associated  with  a  very  flexible  aircraft 
in  flight.  Various  integration  methods  exist  for  linear  structural  dynamics  problems.  However,  a  review 
of  the  literature  indicates  little  material  associated  with  the  integration  of  nonlinear  structural  EOM  of 
relatively  large  order.  Moreover,  for  the  problem  of  simulation  of  very  flexible  aircraft,  a  combination  of 
flight  dynamics  and  aeroelastic  degrees  of  freedom  must  be  integrated  concurrently.  A  modified  first  and 
second  order  Generalized-a  Method  along  with  an  implicit  sub-iteration  scheme  were  developed.  It  has 
shown  good  agreement  with  predictor/corrector  integration  schemes  for  a  reduced  set  of  linear  EOM.  The 
method  is  also  seen  to  be  numerically  stable  when  compared  to  non-dissipative  time  marching  integration 
schemes  and  requires  less  computational  time  compared  to  predictor/corrector  methods  for  the  full  set  of 
nonlinear  EOM. 


II.  Introduction 

Recent  advances  in  airborne  sensors  and  communication  packages  have  brought  the  need  for  high-altitude 
long-endurance  (HALE)  aircraft.  These  platforms  can  be  categorized  under  three  broad  missions  supporting 
either  the  military  or  civilian  community.  The  missions  include  airborne  Intelligence,  Surveillance,  and  Re¬ 
connaissance  (ISR)  for  the  military,1  network  communication  nodes  for  the  military  and  civilian  community,2 
and  general  atmospheric  research.2  Due  to  the  mission  requirements,  the  desired  vehicles  are  characterized 
by  high  aspect  ratio  wings  and  slender  fuselages.  Example  of  mission  optimization  studies  for  this  class  of 
vehicle  can  be  found  in  Ref.  1  where  the  authors  show  that  the  HALE  aircraft  are  required  to  have  a  fuel 
fraction  greater  than  66%,  resulting  in  a  very  small  structural  weight  fraction.  Therefore,  the  combination 
of  high  aerodynamic  efficiency  and  low  structural  weight  fraction  results  in  inherently  very  flexible  vehicles. 
The  HALE  vehicle  may  then  present  large  dynamic  wing  deformations  at  low  frequencies,  presenting  a  direct 
impact  into  the  flight  dynamic  characteristics  of  the  vehicle. 

In  the  process  of  developing  and  subsequently  integrating  the  resulting  set  of  nonlinear  second-  and  first- 
order  differential  equations,3  an  adaptation  of  the  Generalized-a  Method4'5  was  developed  by  the  authors. 
The  resulting  set  of  second-order  nonlinear  elastic  equations  of  motion  (EOM)  and  the  first-order  body 
EOM  are  coupled  with  Peters6’7  finite  state  inflow  model.  The  resulting  set  of  second-  and  first-order 
differential  equations  are  then  integrated  using  a  modified  implicit  Generalized-a  Method.  The  Generalized- 
a  Method  is  a  time  marching  high-frequency  dissipative  integration  scheme  developed  for  linear  systems. 
When  integrating  structural  dynamical  problems,  frequently  the  need  arises  for  a  dissipative  method  to 
prevent  high  frequency  numerical  errors  from  accumulating  and  affecting  the  low  frequency  dynamics  of 
interest. 

The  high-frequency  errors  are  due  to  the  integration  of  a  set  of  stiff  set  of  equations.8,9  Stiff  systems  are 
defined  as  one  with  a  large  condition  number  (ratio  of  the  largest  singular  value  divided  by  the  smallest) 
or  a  system  with  a  very  wide  spread  of  time  constants.10  The  resulting  set  of  differential  equations3  are 
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inherently  stiff.  Additionally  the  equations  can  be  augmented  with  constraint  algebraic  equations  to  impose 
relative  or  absolute  motion  constraints. 

A.  Previous  Work 

Several  authors  have  expanded  the  field  of  structural  dynamics  time  marching  integration  schemes.  New- 
mark11  was  one  of  the  earliest  researchers  who  saw  the  need  for  dissipative  numerical  integration  schemes  for 
structural  dynamics.  That  work  was  followed,  among  others,  by  Hilber,  Hughes,  and  Taylor12  (HHT)  who 
extended  the  Newmark  method.  Several  other  authors  contributed  to  the  body  of  knowledge  are  summarized 
by  Fung’s  review  of  time  marching  algorithms  using  numerical  dissipation.13 

Several  of  the  time  marching  dissipative  methods  were  then  brought  together  in  a  single  second  order 
formulation  by  Chung  and  Hulbert.4  This  formulation  incorporates  the  Newmark,  HHT,  Wilson-0,  and 
trapezodial  methods.  Seeing  the  need  for  a  first-order  high-frequency  dissipative  method,  Jansen,  Whit¬ 
ing,  and  Hulbert5  extended  the  second-order  method  to  a  first-order  one  originally  for  the  integration  of 
computational  fluid  dynamics  equations.  Supporting  the  use  of  numerically  dissipative  integration  schemes, 
Cardona  and  Geradin8  proved  the  importance  of  using  these  schemes  for  the  integration  of  constrained  EOM 
with  finite  rotations. 

Alternate  methods  of  dealing  with  numerical  instability  associated  with  stiff  system  of  equations  have 
been  developed  by  several  researchers.9, 14-17  These  researchers  have  developed  various  momentum  and 
energy  preserving  schemes  as  well  as  momentum  preserving  and  energy  decaying  schemes.  These  methods 
typically  have  slightly  better  convergence  properties  than  the  Generalized-a  Method.  All  the  methods  have 
been  shown  to  work  well  with  conservative  or  state  independent  generalized  forces.  Zhou  and  Tamma18  have 
developed  a  more  general  single  step  integration  scheme  for  linear  structural  dynamics  problems  incorporating 
numerical  dissipation.  Kane  and  Levison19,20  and  Sochet21  extended  the  integration  schemes  by  developing 
various  checking  functions  used  to  post  process  a  numerical  integration  and  determine  its  error. 

Panda22  presented  the  highlights  of  a  Newton-Raphson  sub-iteration  technique  which  utilized  the  parent 
Newmark-/?  time  integration  scheme  for  second-order  EOM  developed  for  flexible  rotorcraft  problems.  The 
main  difference  between  Panda’s22  formulation  and  the  method  presented  here  is  the  incorporation  of  first- 
order  differential  equations.  Additional  differences  with  respect  to  his  paper  are  the  use  of  the  Generalized-a 
time  marching  integration  scheme,  a  generic  beam  model  formulation,  detailed  implementation  schemes,  and 
comparative  results. 

Researchers  studying  nonlinear  aeroelastic  effects  have  used  a  variety  of  techniques.  Patil  and  Hodges  2005 
IFASD,  utilized  a  second-order,  central-difference,  time  marching  algorithm  with  high  frequency  damping, 
page  14.  Drela23  uses  a  second  order  backward  difference  method,  with  a  sub-iteration  step  in  his  ASWING 
code.  Tang  and  Dowell24,25  used  a  reduced  order  model  and  a  Runge-Kutta  integration  scheme.  Patil,  et. 
al.26,27  utilized  a  time-marching  scheme  based  upon  space-time  finite  elements.  Brown28  and  Cesnik  and 
Brown29  utilized  a  trapezoidal  integration  scheme  for  flexible  and  limited  rigid  body  motions.  These  various 
methods  have  shown  difficulty  in  integrating  the  EOM  developed  by  Shearer  and  Cesnik3  either  through 
numerical  instability  or  severely  increased  computational  burden. 

The  Generalized-a  Method  was  selected  based  upon  its  relative  ease  of  implementation  with  the  current 
EOM  modeling  and  the  availability  of  both  first-  and  second-order  formulations.4,5  The  two  methods  are 
modified  using  an  implicit  integration  scheme  detailed  by  Geradin  and  Rixen30  for  nonlinear  second-order 
EOM.  The  second-order  Generalized-a  Method  is  used  to  integrate  the  flexible  EOM,  while  the  first-order 
method  is  used  to  integrate  the  body  EOM  and  remaining  differential  equations  developed  by  Shearer  and 
Cesnik.3  The  use  of  second-  and  first-order  integration  schemes  keep  the  size  of  the  resulting  sub-iteration 
tangent  matrices  significantly  smaller  than  if  the  second-order  equations  were  transformed  to  a  set  of  first- 
order  differential  equations. 

B.  Objective  of  the  paper 

The  objective  of  this  paper  is  to  present  an  implicit  time  marching  numerical  integration  method  for  use  with 
coupled  first-  and  second-order  nonlinear  differential  equations  of  motion,  termed  the  Modified  Generalized-a 
Method.  The  proposed  method  addresses  long  term  integration  stability  and  computational  performance  for 
a  large  nonlinear  elastic  system. 


2  of  20 


American  Institute  of  Aeronautics  and  Astronautics 


III.  Theoretical  Development 


The  theoretical  development  is  comprised  of  four  sub-sections.  The  first  section  is  a  review  of  the 
Generalized-a  Method  for  first-  and  second-order  linear  systems.  The  second  section  presents  a  summary  of 
the  particular  first-  and  second-order  EOM  to  be  solved.  The  third  section  reviews  and  extends  Geradin  and 
Rixen’s  30  method  for  nonlinear  systems  using  a  dissipative  time  marching  integration  scheme  and  provides 
the  details  for  the  very  flexible  aircraft  EOM.  The  final  section  presents  the  details  of  the  convergence  criteria 
required  for  each  time  step. 


A.  Review  of  the  Generalized-a  Method 

The  Generalized-a  Method4  is  designed  to  solve  the  second  order  linear  differential  equation  of  the  form 

Ma  +  Cv  +  Kd  =  F  (1) 


where  M,  C,  and  K  are  generalized  mass,  damping  and  stiffness  matrices,  a,  v,  and  d  are  generalized 
acceleration,  velocity,  and  displacement,  and  F  is  the  generalized  force  vector.  The  Generalized-a  Method 
then  solves  for  the  a  discrete  time  step,  n,  using 


dn+i 

v„+i 

F  (tn+ l-a/2) 

where  h  is  the  time  step  defined  by 
and 


dn  +  hvn  +  h2  ^  Q  -  /?2  j 

vn+h  ((1  -  72)  an  +  72a„+i) 
M&n+l-a^  -f  Gvn+ 1— ar/2  T  Kdn+l—cif2 

h  —  tn+l 


dn+l-a/2 
Vn+1— a/2 

an+l-a,„2 

^n+l-a/2 


(1  -  a/2)dn+i  +a/2dn 

(l-OfaWn+l  +a/2Vn 
(1  “  ttmj)  an+l  +  Qmz&n 
(1  Q/2  )  tn+ 1  4-  ay2tn 


(2) 

(3) 

(4) 


(5) 


(6) 

(7) 

(8) 
(9) 


The  parameters  a/2 ,  am2,  72,  and  /32  are  used  to  control  the  amplification  of  high  frequency  numerical  modes 
which  are  not  of  interest.  If  the  parameters  are  chosen  correctly,  HHT,  Newmark  or  WBZ  methods  can  be 
recovered.  However,  for  this  study  the  following  relationships  are  used 


72 

A 


ah 


1 

=  2  “  am2  +  a/2 
1  2 

^  (1  —  ^7712  ) 

^PoQ2  1 
P0O2  “J”  ^ 

—  PoQ  2 

P0O2  ^ 


(10) 

(11) 

(12) 

(13) 


where  the  subscript  2  refers  to  the  second-order  system,  Eq.  1,  and  will  be  necessary  to  distinguish  72, 
am2,  and  a/2  values  from  the  first-order  EOM  counterpart.  The  single  parameter  is  used  to  control 
the  numerical  dissipation  above  the  normalized  frequency  h/T,  where  T  is  the  period  associated  with  the 
highest  frequency  of  interest  and 

0  <  Poo  <  1  (14) 

If  poo  is  chosen  to  be  unity  then  the  trapezoidal  method  is  recovered.  If  p^  is  chosen  to  be  0  then  frequencies 
above  h/T  will  be  dissipated  in  one  time  step,  a  so  called  “asymptotic  annihilation.” 

In  a  similar  manner  Jansen,  Whiting  and  Hulbert5  developed  the  first-order  Generalized-a  Method  for  a 
first-order  system  of  the  form 

x  =  ax  (15) 
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where  the  following  relationships  are  employed 


»^n+ami  —  ^^n+a/j 

(16) 

xn+\  =  xn  +  hxn  +  h-yi  (in+i  -  xn) 

(17) 

%n+ccmi  =  T  (in^-i  itn) 

(18) 

•Tn+a/j  =  T  Oy,  (-Tn+1  “  -Tn) 

(19) 

In  a  similar  manner  to  the  second-order  Generalized-a  Method  the  free  parameters  71 , 
chosen  in  terms  of  a  single  high  frequency  spectral  radius  parameter,  p ^ ,  by 

ami,  and  aj1  can  be 

1 

7i  g  Qmi  —  aB 

(20) 

1  —  pooi  \ 

mi  2  \  1  +  Pool  ) 

(21) 

Q/l  1  +  Pa O! 

(22) 

B.  Differential  Equations  for  a  Very  Flexible  Aircraft 

The  nonlinear  equations  of  motion  for  a  very  flexible  aircraft  have  been  presented  in  Shearer  and  Cesnik.3 
An  orthogonal  reference  frame  B  is  placed  at  point  O,  which  in  general  is  not  the  center  of  mass  of  the 
vehicle,  as  shown  in  Figure  1 .  The  resulting  set  of  differential  equations  can  be  summarized  as 


Mpp'i  —  —MpsP  —  Cffc  —  CfbP  —  Kppe  +  Rp  (23) 

Mbb/3  =  -Mbf'(  —  CbbP  —  Cbf(  +  Rb  (24) 

C  =  (25) 

Pb  =  [  CGB  0  ]  P  (26) 

A  =  Fl{  l  }+F2{  0  }+FsA  (27) 


where  Eq.  23  is  the  governing  nonlinear  structural  second  order  EOM,  Eq.  24  is  the  B  reference  frame  first 
order  nonlinear  EOM,  Eq.  25  is  the  propagation  of  the  orientation  of  the  B  reference  frame  using  quaternion 
parameters,10  Eqs.  26  and  27  represents  the  unsteady  aerodynamic  effects  through  induced  inflow  over  the 
lifting  surfaces.  The  variables  are  defined  as 

e  =  strain  vector 

P  =  vector  of  translational  and  rotational  velocities 
C  =  quaternion  parameters  for  B  reference  frame  orientation 
Pb  =  vector  components  of  B  reference  frame  location 
A  =  unsteady  inflow  velocities 

Three  possible  solutions  to  Eqs.  23-27  are  given  in  Shearer  and  Cesnik.3  First,  by  reducing  the  order  of 
the  equations  by  eliminating  all  elastic  DOF,  a  simple  rigid  body  solution  emerged.  Second  is  a  linearized 
solution,  where  the  generalized  mass  matrix  is  only  a  function  of  the  initial  state.  Third  is  a  full  nonlinear 
simulation  where  the  generalized  mass  matrix  is  updated  at  each  time  step.  The  matrices  of  Eqs.  23-27  are 
functions  of  the  states,  e,  (3,  £,  A  as 

Mff,  Mfb,Mbf,Mbb 
Cff,Cbf 
Cfb,Cbb 

C°B 
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Inertial  Ftim«  (G) 


Figure  1.  Basic  Body  Reference  Frame  and  Vehicle  Coordinates 


C.  Implicit  Integration  Scheme  Utilizing  Generalized-a  Method 

The  Generalized-a  Method  was  designed  for  linear  systems  but  is  implemented  here  in  solving  a  non-linear 
problem.  In  general  all  numerical  integration  schemes  follow  the  flow  given  in  Figure  2.  The  differences 
between  integration  schemes  are  imbedded  in  the  “Subiteration  Routine”  block.  For  a  trapezodial  method, 
this  block  consists  of  the  amplification  matrix,  A,  which  solves 


xn+1  =  [A]  xn  (29) 

where  n  is  a  discrete  time  step  and  x  are  the  system  states.  In  the  modified  Generalized-a  Method  the 
basic  concept  at  each  time  step  is  to  predict  the  states  and  their  time  derivatives  and  employ  a  sub-iteration 
Newton-Raphson  method  to  correct  the  state  predictions.  The  sub-iteration  is  repeated  until  a  user  defined 
tolerance  is  met.  The  implicit  integration  scheme  chosen  resembles  Geradin’s  and  Rixen’s  method30  and  the 
flow  is  shown  in  Figure  ;i.  The  specific  convergence  flow  is  shown  in  Figure  4. 


1.  Predictors 

To  begin  the  sub-iteration  loop,  Figure  3,  the  states  at  time  step  n  +  1  are  predicted.  For  the  second  order 
EOM,  Geradin  and  Rixen30  provide  a  set  of  predictors  <j*+1 ,  <j*+1 ,  and  q*+1  given  the  states  at  the  current 
time  step,  n,  qn,  qn,  and  qn  as 

9n+ 1 

<7n+l 
9n+l 

In  a  similar  manner  to  Geradin  and  Rixen,30  the  first  order  EOM  predictors,  x*+)  and  Xn+i  are  proposed 
as 

^n+l  =  't'n  T  h  (1  7l )  Xn  (33) 

K+i  =  0  (34) 


—  qn  +  hqn  + 


6-' 


#2  h2qn 


qn  +  (1  ~  72)  hq„ 

0 


(30) 

(31) 

(32) 
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Figure  2.  Basic  Numerical  Integration  Flow 
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Figure  3.  Sub-iteration  Routine 
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2.  Newton- Raphson:  Residual  Terms  and  the  Tangent  Matrix 

The  predictors  of  Eqs.  30-34  are  substituted  back  into  the  governing  differential  equations,  Eqs  23-27,  where 
all  the  terms  are  moved  to  the  left  side.  The  result  is  a  set  of  residual  terms,  r,  defined  as 


rf  =  Mppe  +  Mfb/3  +  Cff£  +  CfbP  +  Kffc  —  Rf  (35) 

r;,  =  MbbP  +  Mbf£  +  CbbP  +  Cbfz  —  Rb  (36) 

rq  =  C  +  ^cC  (37) 

rP  =  Pb  —  [  CGB  0]p  (38) 

n  =  A-Fi  |  |  j-F2  j  j  j-F3A  (39) 

where 

r  =  [  r/  rj  r,  rp  r,  j  (40) 

A  Taylor  series  expansion  of  the  residual  vector,  ,  at  time  step  n  +  1  and  sub-iteration  step  k  +  1  yields 

rftl  =  rkn+ 1  +  []fqrkn+i]  (fl&l  ~  <£+i)  +  W.O.T.  (41) 


and 

9  =  [  eT  /?T  CT  VTb  At]T  (42) 

Setting  the  higher  order  terms  to  zero,  assuming  r*%\  =  0,  and  defining  the  tangent  matrix  S^+1 
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where 


A  qk=qknl\-qkn+i 


(44) 


A qk  is  solved  as 

A  =  (45) 

Here  A qk  can  be  the  correction  term  for  either  a  first  or  second  order  differential  equation  depending  upon 
the  structure  of  the  tangent  matrix.  For  the  governing  differential  equations,  Eqs.  23-27,  the  full  tangent 
matrix  takes  the  form 

[eSF(r/)n+l]  [a3^(r/)n+l]  [a^(r/)n+l]  [dApJ,  (r/)n+l]  (r/)n+l] 

[aAeF(r’6)n+l]  [s^(rf>)n+l]  [aS^Wn+l]  [a/fp^  (ri>)n+l]  ['&£& (r&)£+l] 

[a^^(r«)n+l]  [a^'(r9)n+l]  [d<\pkD  (r9)n+l]  [sAA^  (r<l)n+l] 

[aAeF(rp)n+l]  [flA/jt  (rp)n+l]  [s3?I’(rl,)«+l]  [a/fp*  (rp)n+l]  [aAA^  (rp)n+l] 

[s^(r*)n+l]  [sA/3l  (r<)n+l]  \_dE<*  (ri)n+l]  [a^fp*  (ri)n+l]  UaA*-'  (r«)n+l] 

The  tangent  matrix  S„+1  given  by  Eq.  40  in  practice  tends  to  have  a  large  condition  number  due  to 
the  stiffness  of  the  governing  differential  equations.  This  creates  a  problem  with  numerical  accuracy  when 
inverting  S*+1 .  In  order  to  reduce  the  condition  number  and  improve  numerical  accuracy,  the  tangent  matrix 
is  scaled  with  the  scalar  quantities  dj  according  to 

Skn+Uj=Skn+Uj.dj  (47) 

The  scalar  values,  dj,  are  found  so  to  minimize  the  condition  number  of  S*+1.  Equation  45  is  then  modified 

as 

A<f  =  -[S*+I]_1r*+1  (48) 

A Qk  =  [  diAgf  ofeAgf  •••  dsAc/f  ] 

3.  Correction  Terms 

For  second  order  differential  equations  the  states  and  their  derivatives  are  updated  as30 

QnX\  =  ?n+l+A?fc 

tit  1  =  ti+i  +  (^)  (49) 

till  =  + 

For  the  first  order  differential,  equations  the  updated  states  and  the  first  derivative  terms  are  found  to  be 

till  =  ti+i  +  &ti 

till  =  ti+i  +  (~)*ti  (50) 

This  is  derived  starting  with  the  update  equation  of  qn+\ ,  Ref.  5,  Eq.  13 

9n+i  =  Qn  +  hqn  +  hji  (<7n+i  -  qn)  (51) 
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where  the  following  substitutions  were  made 

q  =  y,  h  =  At,  71  =  7 

Equation  51  can  be  solved  for  qn+\  as 


Substituting  Eq.  33  into  Eq.  53 


Qn+ 1  =  (qn+i  ~  {qn  +  h  (1  -  71)  qn)) 


9n+l  —  ^  (Qn+i  <7n+l) 


and  defining 

A Q  =  Qn+ 1  —  Qn+l  (55) 

Eq.  50  is  recovered.  Using  Eqs.  49  and  50  the  partial  derivatives  of  the  states  and  their  derivatives  with 
respect  to  A q  can  now  be  computed.  For  second  order  differential  equations  the  derivatives  are 


'  d  .; 

1 

9  A  qq 

'  9  ; 

72 

9  A  qq_ 

fah 

'  9  ‘ 

=  1 

_dAqq 

and  for  first  order 


—  -1  =  1 
9  A  qq  71  h 

d  1  _ 
dAqq 


Using  Eqs.  23,  46,  56  and  58  the  first  row  of  the  tangent  matrix  is  given  as 


(‘S'n+^First  Row 


[(0^*)  Mff  +  (^s)  Cff  +  Rff] 
[(^)  Mfb  +  Cfb  _  [e^]] 
[&**■] 


where  the  matrices  Mff,Cff ,  Mfb,  and  Cfb  are  assumed  to  be  constant  for  the  derivation  of  the  tangent 
matrix. 

D.  Convergence  Criteria 

Figure  4  presents  the  high  level  flow  of  the  convergence  routine.  The  Newton-Raphson  method  outlined 
above  works  well  for  the  majority  of  time  steps  and  sub-iteration  steps.  However,  given  the  large  number 
of  states  that  are  being  solved,  the  Newton-Raphson  on  occasion  will  not  yield  a  lower  norm  of  the  residual 
vector.  If  that  happens  the  state  correction,  A q  is  modified  using  a  line  search  algorithm,  such  that 


Aq  =  a(-[S*+1]  \kn+1) 
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Crude  bracket  values  for  the  scaling  parameter  a  are  first  found  by  calculating  the  residual  vector  at  various 
values  of  a  such  that 

0<Q(<a<au<l  (62) 

The  lower  and  upper  bounds  (subscripts  l  and  u)  on  a  are  then  used  in  a  quadratic  curve  formula7 


(  {<*m  ~  <*l)  INI  +  {al  ~  al)  1km |1  +  (of  ~  Qm)  \K\\\ 

\  («*rn  -  o„)  llnll  +  (a«  -  at)  ||rm||  +  (a;  -  am)  ||r„||  J 


(63) 


where 


Cu  —  Oil 
2 


(64) 


and  rm  is  evaluated  at  am.  Equation  63  is  iterated  upon  until  a  satisfactory  convergence  on  a  is  reached. 
The  change  in  the,  A q  is  then  updated  in  accordance  with  Eqs.  49,  50,  and  61.  In  almost  all  cases  this  line 
search  method  provides  an  excellent  update  state,  A q. 

It  is  possible  however  for  the  Newton-Raphson  method  to  converge  to  a  local  minimum  of  ||r||2-  In  this 
case,  the  line  search  scaling  parameter,  a  will  be  zero.  This  situation  is  determined  by  monitoring  the  value 
of  ct  for  several  sub-iteration  steps.  If  a  <  e  for  more  than  a  user  defined  number  of  sequential  sub-iteration 
steps,  than  a  is  arbitrarily  set  to 

a  =  0.25  +  0.25 h  (65) 


where  fc,  is  the  number  of  times  a  local  minimum  has  been  reached  and  £  is  a  user  defined  and  «  0.  The 
state  update,  Aq  is  then  computed  using  Eq.  61 .  Using  this  heuristic  approach,  the  Newton-Raphson  search 
is  moved  away  from  a  local  minimum  and  allowed  to  continue  searching  for  the  global  minimum  of  ||r||2.  In 
practice  this  method  has  shown  excellent  results  at  resolving  convergence  to  local  minimums. 


IV.  Numerical  Examples 


Two  different  models  are  presented  here  to  highlight  the  main  characteristics  of  the  proposed  method. 
The  first  model  is  a  simple  cantilevered  beam  shown  in  Figure  5  with  properties  given  in  Table  I . 


t 


/ 
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/ 
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Figure  5.  Cantilevered  Beam 

Three  integration  techniques  were  run  on  the  beam  model:  Matlab’s  ODE15S,  Trapezoidal  Method, 
Modified  Generalized-a  Method  .  For  the  beam  model  Eq.  23  is  solved  where  the  B  reference  frame  states, 
(3,  are  removed.  Linear  and  nonlinear  solutions  are  presented  by  either  setting  Mff  and  Cfj  to' the  initial 
state  or  being  updated  at  each  time  step.  A  plot  of  the  linear  solution  with  the  three  integration  methods 
is  shown  in  Figure  6  and  a  table  of  the  relevant  features  is  presented  in  Table  2.  All  three  methods  are 
seen  to  provide  similar  results.  The  trapezoidal  method  and  Generalized-a  Method  both  provide  a  relatively 
quick  solution  to  the  differential  equations,  while  Matlab’s  ODE15S  takes  almost  2  orders  of  magnitude 
longer  due  to  the  stiffness  of  the  equations.  For  the  nonlinear  case,  Figure  7,  Matlab’s  ODE15S  is  seen  to 
diverge  from  the  other  solutions.  This  is  due  to  the  inability  of  the  Matlab’s  solvers  to  handle  high  accuracy 
stiff  ODEs.  Tighter  tolerances  did  move  Matlab’s  ODE15S  solution  slightly  closer  to  the  trapezoidal  and 
Modified  Generalized-a  Method  solutions,  but  at  the  risk  of  early  termination  or  increased  integration  time 
of  several  orders  of  magnitude.  Matlab’s  other  solvers  performed  in  a  similar  manner.  This  is  not  surprising 
as  Matlab’s  ODE  solvers  are  intended  for  non-stiff,  lower  order  ODEs.  Table  3  provides  a  comparison  of 
computation  time  for  the  various  integration  methods  and  cases. 
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Cantilevered  Beam  Properties 

Property 

Value 

Units 

Beam  Length 

1.0 

m 

Kn  =  EA 

106 

N 

K22  =  GJ 

50 

N  •  m 2 

A'33  =  EI2 

50 

N-m2 

Ku  =  EI3 

103 

N-m 2 

Mass  per  unit  length 

0.2 

kg  ■  m_1 

In  per  unit  length 

10“4 

kg  ■  m 

I22  per  unit  length 

10“6 

kg  ■  m 

I33  per  unit  length 

10~4 

kg  ■  m 

Elements  in  beam 

10 

— 

Number  of  Second-Order  states 

40 

— 

Table  1.  Geometric,  Stiffness,  and  Inertia  Properties  of  Cantilevered  Beam  (only  non-zero  terms  are  listed) 


Figure  6.  Linear  Elastic  Solution  of  Cantilevered  Beam  with  Tip  Force  Actuation  of  10sin(20t) 


Linear  Elastic  Solution  of  Cantilevered  Beam,  Numerical  Results 

Integration  Scheme 

CPU  Time 

Maximum  Difference 

Convergence  Parameters 

MatlabODE15S 

1746.80  s 

0 

Relative  =  10~3,  Absolute  =  10~5 

Trapezoidal 

23.64  s 

0.0064 

NA 

Generalized-a  Method 

28.47  s 

0.0124 

IISIBOliUeBBHI 

Table  2.  Comparative  Results  of  Various  Integration  Schemes  for  a  Cantilevered  Beam 
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Figure  7.  Nonlinear  Elastic  Solution  of  Cantilevered  Beam  with  Tip  Force  Actuation  of  10sin(20t) 


Computational  Times  for  Cantilevered  Beam 

Integration  Scheme 

Linear  Solution 

Nonlinear  Solution 

MatlabODE15S 

1746.812s 

9946.000s 

Trapezoidal 

23.641s 

38.781s 

Modified  Generalized-Q  Method 

28.469s 

80.032s 

Table  3.  Computational  Times  for  Trapezoidal  and  Modified  Generalized-a  Method  Integration  Schemes  for 
Very  Flexible  Aircraft 
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For  the  next  result,  the  cantilevered  beam  is  pinned  at  the  50%  location  as  seen  in  Figure  8.  In  this  case 


z 

t 


Figure  8.  Cantilevered  Beam  Constrained  at  Mid  Point 


the  elastic  equations  of  motion,  Eq.  23,  have  been  modified  by  Su  and  Cesnik31  and  the  authors  through 
the  stiffness  matrix  to  include  algebraic  constraint  equations  with  Lagrange  Multipliers,  A £.  The  augmented 
system  given  in  Eq.  23  without  the  B  reference  frame  linear  and  angular  velocities,  /?, becomes 


where  $,  contains  the  constraint  relation,  Xi  is  the  Lagrange  multiplier,  *  is  a  time  index,  and  hi  and  ho 
are  the  displacement  information  at  the  current  time  step  i  and  the  initial  time  step  0.  For  this  case  Matlab 
does  not  have  a  solver  which  can  handle  differential  algebraic  equations  (DAEs)  higher  than  index  1.  Recall 
the  index  of  a  DAE  is  defined  as  the  number  of  times  which  the  algebraic  equation  must  be  differentiated 
before  a  standard  ODE  form  is  reached.  For  the  current  case  the  constraint  equation  (before  expressing  it 
in  discrete  form)  is  of  the  form 

/(£)  =  0  (67) 

Differentiating  /(e)  two  times  yields 


such  that  the  later  is  in  ODE  form.  The  results  for  this  case  with  a  linear  solution  are  seen  in  Figure  9. 
While  the  trapezoidal  case  appears  to  be  stable  for  the  linearized  solution,  a  closer  examination  of  some  of 
the  discrete  eigenvalues  of  the  amplification  matrix,  Table  4,  reveal  a  slight  instability  (greater  than  unity 
or  repeated  eigenvalues  on  the  unit  circle) .  The  eigenvalues  with  only  real  parts  and  associated  eigenvectors 
are  due  to  the  Lagrange  multipliers.  By  examining  the  unstable  elastic  eigenvalues  it  is  found  that  they  are 
also  controlled  by  the  Lagrange  multipliers.  The  slight  instability  is  not  seen  over  relatively  short  periods 
of  integration.  This  is  consistent  with  the  proof  of  Cardona  and  Geradin.8  For  the  nonlinear  solution, 
Figure  10,  it  is  seen  that  the  trapezoidal  method  is  unstable.  This  can  also  be  seen  from  the  residual  term  as 
shown  in  Figure  If.  Also  from  Figure  10,  the  Generalized-a  Method  maintains  long  term  stability.  However 
due  to  the  high  frequency  numerical  damping  of  the  Modified  Generalized-a  Method,  there  is  a  loss  of  high 
frequency  content.  Solution  time  is  also  longer  for  the  Modified  Generalized-a  Method  ,  Table  5,  due  to  the 
recursive  nature  of  the  sub-iteration  scheme  until  a  satisfactory  residual  term  is  obtained.  The  user  must 
trade  long  term  stability  with  the  loss  of  high  frequency  content. 

For  very  flexible  aircraft  repeated  eigenvalues  on  the  unit  circle  can  come  from  the  aircraft  configuration 
(joined  wing  concept),1  unconstrained  rigid  body  degrees  of  freedom  due  to  a  free  flying  aircraft,  or  an  aircraft 
controller.  A  model  based  upon  Ref.  3  is  shown  in  Figure  1  and  relevant  physical  properties  are  summarized 
in  Table  6.  Here  results  are  presented  for  a  nonlinear  flexible  simulation  where  the  aircraft  is  given  a  square 
aileron  input,  Figure  12.  Representative  solutions  of  the  longitudinal  and  vertical  B  reference  frame  velocities 
are  shown  in  Figure  13  and  pitch  and  yaw  rates  in  Figure  14  for  a  trapezoidal  and  Generalized-a  Method 
integration.  While  the  two  methods  track  resonably  well  for  the  first  10  seconds  and  then  the  trapezoidal 
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time,  s 


Figure  9.  Linear  Elastic  Solution  of  Pinned  Cantilevered  Beam  with  Tip  Force  Actuation  of  10sin(20t) 


Eigenvalues  due  to  Lagrange  Multiplier  Constraints 

Eigenvalue 

Real  Part 

Imaginary  Part 

State 

^75 

-0.999992847 

0 

Lagrange  Multiplier 

•^78 

-0.999996417 

0 

Lagrange  Multiplier 

•^81 

-0.999999461 

0 

Lagrange  Multiplier 

•^76,77 

-1.000003576 

+/-  0.000006194i 

Strain 

•^79,80 

-1.000001791 

+/-  0.000003103i 

Strain 

ihim 

-1.000000269 

+/-  0.000000467i 

Strain 

Table  4.  Comparative  Results  of  Various  Integration  Schemes  for  a  Cantilevered  Beam 


Figure  10. 


Nonlinear  Elastic  Solution  of  Pinned  Cantilevered  Beam  with  Tip  Force  Actuation  of  10sin(20£) 
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Figure  11.  Residual  Term,  ||r||2,  of  Pinned  Cantilevered  Beam  with  Tip  Force  Actuation  of  10sin(20<) 


Computational  Times  for  Cantilevered  Beam 

Integration  Scheme 

Linear  Solution 

Nonlinear  Solution 

Trapezoidal 

73.875s 

110.531s 

Modified  Generalized-a  Method 

75.156s 

258.375s 

Table  5.  Computational  Times  for  Trapezoidal  and  Modified  Generalized-a  Method  Integration  Schemes  for 
Very  Flexible  Aircraft 


Figure  12.  Open  Loop  Aileron  and  Rudder  Commands 
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Model  Parameters 

Property 

Value 

Units 

Light 

Heavy 

Fuselage  Length 

26.4 

m 

Wing  Span 

58.6 

m 

Wing  Area 

196.3 

m2 

Root  Chord 

4.5 

m 

Tip  Chord 

2.2 

m 

Aspect  Ratio 

17.5 

— 

Horizontal  Tail  Span 

18.0 

m 

Horizontal  Root  Chord 

3.5 

m 

Horizontal  Tip  Chord 

2.45 

m 

Vertical  Tail  Span 

4.0 

m 

Vertical  Root  Chord 

2.45 

m 

Vertical  Tip  Chord 

2.0 

m 

Wing/Horizontal  Tail  Airfoil 

NACA  4415 

— 

Vertical  Tail  Airfoil 

NACA  0012 

— 

Aileron  Location 

16.3  to  22.8 

m 

Aileron,  Elevator,  Rudder  Chord 

0.2  CioCai 

Elevator  Location 

1.8  to  9.0 

m 

Rudder  Location 

0.8  to  3.2 

m 

Aircraft  Angle  of  Attack 

0.64° 

6.37° 

— 

Elevator  Deflection  Angle 

-4.11° 

-13.43° 

— 

Fuel  Mass 

0 

32,000 

kg 

Total  Mass 

1.52  •  104 

4.72  •  104 

kg 

Fuel  Fraction 

0.0 

67.8 

% 

TSS  * 

1XX 

9.61  •  105 

1.17 -106 

kg  •  m2 

TSS 

w 

8.21  •  105 

2.94  •  106 

kg  •  m2 

TSS 

1 ZZ 

1.75  •  106 

3.93  ■  106 

kg  •  m2 

0 

0 

kg  •  m2 

TSS 

1XZ 

0 

0 

kg  •  m2 

-1.65- 104 

-4.72  •  104 

kg  •  m2 

Elements  per  wing 

9 

— 

Elements  per  horizontal  tail 

5 

— 

Elements  per  vertical  tail 

5 

— 

Elements  in  fuselage 

10 

— 

Total  Number  of  Elements 

48 

— 

*Note:  Iss  are  the  inertia  properties  in  a  steady  state  configuration 


Table  6.  Geometric  and  Inertia  Properties  of  the  Flexible  Aircraft  Model 
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Pitch  Rate  Longitudinal  Velocity,  m/s 


Figure  13.  B  Reference  Frame  Velocities  for  Trapezoidal  and  Generalized-o  Method  Integration 


Figure  14.  B  Reference  Frame  Velocities  for  Trapezoidal  and  Generalized-o  Method  Integration 
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method  begins  to  diverge.  The  divergence  can  be  better  seen  by  examining  the  two  norm  of  the  residual 
shown  in  Figure  15,  while  the  Modified  Generalized-a  Method  was  commanded  to  keep  the  norm  below 
1.0  (figure  not  shown).  The  penalty  paid  for  this  long  term  stability  is  in  computational  time  as  shown  in 


Figure  15.  Residual  Term,  ||r||2,  of  Trapezoidal  Integration 

Table  7. 


Computational  Times 

Trapezoidal 

3902.8  s 

Modified  Generalized-a  Method 

42261.2  s 

Table  7.  Computational  Times  for  Trapezoidal  and  Modified  Generalized-o  Method  Integration  Schemes  for 
Very  Flexible  Aircraft 


V.  Conclusion 

The  proposed  integration  method,  i.e.,  the  Modified  Generalized-ct  Method,  shows  good  correlation  with 
existing  integration  schemes  for  systems  which  are  stiff  and  have  a  large  number  of  states.  Its  main  limiting 
factors  are  an  increase  in  computational  time  over  simpler  first  order  methods  and  the  attenuation  of  high 
frequency  data  due  to  the  dissipative  nature  of  the  integration  scheme.  It  was  also  shown  that  the  method 
handles  DAE  of  index  higher  than  1  and  preserves  long  term  stability  when  solving  nonlinear  elastic  EOM. 

NASA  Acknowledgement? 
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